This short tutorial will demonstrate some capabilities of the maximus48 Python package on how to:
For some details about the functions used, please refer to the desciption of functions in the corresponding .py files of the package.
The package can be downloaded with: https://pypi.org/project/maximus48/
SSIM is a well-known metrix to find a similarity between two images. It is nicely implemented in skimage: https://scikit-image.org/docs/dev/auto_examples/transform/plot_ssim.html However, in my package I slightly reduced some unnecessary features to make the algorithm work faster and be suitable for parallel processing
Let us start with a simple example with 2D image.
import numpy as np
from maximus48 import SSIM_131119 as SSIM
from maximus48 import var
import tifffile
import matplotlib.pyplot as plt
#set the ROI of image first, the logic corresponds to FIJI (to be read (x,y,x1,y1 at the image - inverse to numpy!)
ROI = (20,20,2048,2048)
folder = '/Users/mpolikarpov/data_server/tomo/'
image_name = 'Platy-Platy_9902_1_01800.tiff'
image2_name = 'Platy-Platy_9902_1_01801.tiff'
#let's read the image and crop it
image = tifffile.imread(folder+image_name)
image = image[ROI[1]:ROI[3], ROI[0]:ROI[2]]
var.show(image)
# please note - we used the built-in function var.show(), it is a bit modified version of matplotlib.pyplot.show()
# this is phase-contrast X-ray image, as you may guess from the contrast
# now let's read all flatfield files from the folder
# in our case flatfield files start with the prefix 'ff_'
#create the list all images in the folder
imlist = var.im_folder(folder)
#read ff-files
flatfield = np.asarray([tifffile.imread(folder+im) for im in imlist if im.startswith('ff')])
flatfield = flatfield[:,ROI[1]:ROI[3], ROI[0]:ROI[2]]
# now it a 3D array of 30 images, stacked along 0-direction
print('the initial shape of flatfield-array is', flatfield.shape)
# each flatfield image looks like this, being however slightly different due to beam instabilities
var.show(flatfield[0])
# please transpose the ff-array for the further ff-correction
flatfield = np.transpose(flatfield, (1,2,0))
print('the final shape of flatfield-array is', flatfield.shape)
# let's divide our image by the flatfield, using SSIM metrix
# images should be set as special classes:
image_class = SSIM.SSIM_const(image)
ff_class = SSIM.SSIM_const(flatfield)
# then, you can calculate SSIM metrics for each pair (data-image) - (ff-image)
index = SSIM.SSIM(image_class, ff_class).ssim()
print('indexes of SSIM are:\n ', index)
# now, simply divide your image by flatfield-image with highest SSIM-index and get a corrected image:
result = image/flatfield[:,:,np.argmax(index)]
var.show(result)
# compare with just a random division by a random ff-image below
bad_result = image/flatfield[:,:,0]
var.show(bad_result)
See, how worse the filtering is (ringy artifacts and horizontal stripes). In addition, the filtration is bad at the edges of the image
Now, let's assume that you have a 3D array of images. The proper way to handle this would be to follow the similar logic as before, however implementing some parallelization. While a simple for loop will be the easiest solution, here I will show the case of using python multiprocessing.
import os
from multiprocessing import Pool
from maximus48.tomo_proc2 import F, tonumpyarray
from contextlib import closing
# fdon't forget to switch off all low-level C-based parallelisation from numpy etc. - gonna be mess otherwise
os.environ['OMP_NUM_THREADS'] ='1'
os.environ['OPENBLAS_NUM_THREADS'] = '1'
os.environ['MKL_NUM_THREADS'] = '1'
We will be using a standard Python high-level parallelisation with Pool. Some things to consider there:
# allocate memory to store flatfield. I do it with a custom class which is based on multiprocessing.Array()
shape_ff = ROI[3]-ROI[1], ROI[2]-ROI[0], 30
ff_shared = F(shape = shape_ff, dtype = 'd')
# read ff-files to memory
ff = tonumpyarray(ff_shared.shared_array_base, ff_shared.shape, ff_shared.dtype)
flatfields = np.asarray([tifffile.imread(folder+im) for im in imlist if im.startswith('ff')])[:,ROI[1]:ROI[3], ROI[0]:ROI[2]].transpose(1,2,0)
for i in range(flatfields.shape[2]):
ff[:,:,i] = flatfields[:,:,i]
# calculate a class of flatfield-related constants (skipped shared memory because it is rather small array)
ff_con = SSIM.SSIM_const(ff)
# initializer to define which variables will be global
def init():
global proj, ff_shared, ff_con, ROI, image_paths
# main function
def Processor(j):
"""
j: int
an index of the file that should be processed
Please note, j always starts from zero
To open correct file, images array uses images[i][j + N_start-1]
"""
#set local variables (unwrap array from shared memory)
ff = tonumpyarray(ff_shared.shared_array_base, ff_shared.shape, ff_shared.dtype)
proj_loc = tonumpyarray(proj.shared_array_base, proj.shape, proj.dtype)
#read image and do ff-retrieval
im = tifffile.imread(image_paths[j])[ROI[1]:ROI[3], ROI[0]:ROI[2]]
index = SSIM.SSIM(SSIM.SSIM_const(im), ff_con).ssim()
#save result to the shared memory
proj_loc[j] = im/ff[:,:,np.argmax(index)]
# set paths of data files
image_paths = (folder+image_name, folder+image2_name)
# allocate memory to store filtered files (only two in my example)
N_files = len(image_paths)
proj = F(shape = (N_files, shape_ff[0], shape_ff[1]), dtype = 'd' )
# Finally, make it all work
cp_count = 3
with closing(Pool(cp_count, initializer = init)) as pool:
pool.map(Processor, np.arange(len(image_paths)))
Let's see what we got
out = tonumpyarray(proj.shared_array_base, proj.shape, proj.dtype)
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 4),
sharex=True, sharey=True)
ax = axes.ravel()
ax[0].imshow(out[0], cmap='gray', vmin=0, vmax=1)
ax[0].set_title('Image 1')
ax[1].imshow(out[1], cmap='gray', vmin=0, vmax=1)
ax[1].set_title('Image 2')
plt.tight_layout()
plt.show()
This is it. You can scale it to big data.
This section is dedicated to simple phase-retrieval with single-distance algorithms. I highly recommend this book to get a comprehensive understanding of maths behind: " T.Salditt et al 2017. Biomedical Imaging. Principles of Radiography, Tomography and Medical Physics "
pixel = 0.1625 * 1e-6 # pixel size
distance = 6.25 * 1e-2 # distances for your measurements
energy = 18 # photon energy in keV
wavelength = var.wavelen(energy)
fresnelN = pixel**2/(wavelength*distance)
# I will use the image from the previous example
filtered = out[0]
import maximus48.sidi_phare as phare
# Paganin algorithm
beta_delta_paganin = 0.1 #0.5
test_paganin = phare.Paganin(filtered, fresnelN, beta_delta_paganin)
var.show(test_paganin)
Note, that Paganin reconstruction retrieves the phase but blurs the image. There is possibility to use unsharp filter to make image sharper (check the formalism here - https://scripts.iucr.org/cgi-bin/paper?mo5009):
w = 2 #width of gauss kernel
stabilizer = 0.3 # weight
proba = phare.anka_filter(test_paganin, w, stabilizer)
var.show(proba)
See, it became sharper because we used high frequencies from the original image
#Modified Bronnikov Algorithm
beta_delta_MBA = 0.1
test_MBA = phare.MBA(filtered, fresnelN, beta_delta_MBA)
var.show(test_MBA)
# Bronnikov Aided Correction
# gamma is a weight of the filter - the same logic as for unsharp filter in Paganin
gamma = 0.1
test_BAC = phare.BAC(filtered, fresnelN, beta_delta_MBA, gamma)
var.show(test_BAC)
# CTF retrieval
from maximus48 import monochromaticCTF as CTF
beta_delta = 0.15
zero_compensation = 0.1 #parameter to cut-off higher frequencies
test_CTF = CTF.single_distance_CTF(filtered, beta_delta, fresnelN, zero_compensation)
var.show(test_CTF)
As you may see, it gives nice results, comparable to Paganin filter in terms of speed and phase retrieval. However, when used with a single distance, this filter introduces stripe artifacts (see below). The best way to minimize them is to use at least 4-distance phase retrieval. Check my other example to get understanding of how it can be done.
var.show(test_CTF[1550:,400:1000])